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Abstract 

We describe and illustrate a simple procedure for identifying a liquid interface from atomic 
coordinates. In particular, a coarse grained density field is constructed, and the interface is defined 
as a constant density surface for this coarse grained field. In applications to a molecular dynamics 
simulation of liquid water, it is shown that this procedure provides instructive and useful pictures 
of liquid-vapor interfaces and of liquid-protein interfaces. 



The Interface 



Definitions of soft-matter interfaces at a molecular level can be ambiguous [H |2] . Due to 
molecular motions, interfacial configurations change with time, and the identity of molecules 
that lie at the interface also change with time. Generally useful procedures for identifying 
interfaces must accommodate these motions. Here, we present a simple and intuitive pro- 
cedure for doing so. The procedure is based upon spatial coarse graining, it applies to 
reasonably arbitrary geometries, and it can be applied at any point in time so that it can 
be used to interpret time dependent phenomena and fluctuations. We find the procedure to 
be useful in a variety of contexts, a few of which are illustrated in this and the next section. 

The basic idea begins with the instantaneous density field at space-time point r, t, 

p{r,t) = J2S{r-r,{t)), (1) 

i 

where rj(t) is the position of the ith particle at time t, and the sum is over all such particles 
of interest. Rendering this field directly provides only vague impressions of interfaces. A 
more manageable field can be formed through coarse graining. Our choice of spatial coarse 
graining is a convolution with the normalized Gaussian functions 

0(r;O = (2vre^)-"/'exp(-rV2e^), (2) 

where r is the magnitude of r, ^ is the coarse graining length, and d stands for dimensionality. 
Applied to p(r, t) we have the coarse grained density field 

p(r,t) = 5^0(|r-r,(t)|;O. (3) 

i 

The choice of ^ will depend upon the physical conditions under considerations. With ^ set, 
we define interfaces to be the {d — l)-dimensional manifold r = s for which 

p(s, t) = c, (4) 

where c is a constant. In other words, we define instantaneous interfaces to be points in 
space where the coarse grained density field has the value c. This coarse grained density 
changes with time as molecular configurations change with time, i.e., s = s{t) = s({rj(t)}). 

For a given molecular configuration, {vi(t)}, Eq. Q can be solved quickly through inter- 
polation on a spatial grid. [T] illustrates what is found for one configuration of a slab of liquid 
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FIG. 1: (a)A snapshot of a slab of liquid water with the instantaneous interface s rendered as a blue 
mesh on the upper and lower phase boundary. The slab is periodically replicated in the horizontal 
directions. (b)The time correlation function governing the spatial fluctuations in the intrinsic 
interface s. Here angle brackets represent an equilibrium average and 5sz{t) = (s(t) ■ z — (s) ■ z) 
where s(t) is the position of the interface at time t and z is the unit vector in the z direction (as 
is indicated in Panel (a)). 

water at conditions of water- vapor coexistence. Details of our simulations are described be- 
low in the Methods section. We have taken {ri(t)} to refer to the positions of all the oxygen 
atoms in the system, and because the bulk correlation length of liquid water is about one 
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FIG. 2: Snapshots of an assembling pair of melittin dimers. (a)Snapshot at t = showing exphcit 
solvent and melittin dimers rendered as green and yellow ribbons, (b) Snapshot at t = excluding 
the explicit solvent with the instantaneous interface s rendered as a blue mesh, (c) Snapshots of 
interface at t = 240ps, t = 377ps, and t = 743ps. Red, blue, and white bands illustrate cross- 
sections in s as seen from the top in (d). 

molecular diameter, we have used use ^ = 2.4A; further, we have used c = 0.016A~'^, which 
is approximately one-half the bulk density. 

The pictured instantaneous water-vapor interface s resembles a wavy sheet dividing the 
liquid and vapor-like regions. In accordance with the predictions of capillary wave theory |3l 
m |5], the undulating height fluctuations in s exhibit transverse correlations that extend 
across the system. As shown in[T], for the size system pictured, these fluctuations de-correlate 
on time-scales of tens of picoseconds. 

This approach belongs to a broader class of interface identification algorithms that build 
upon the assumption that the interface is well described as a molecularly sharp c? — 1 di- 
mensional manifold that is made rough by collective thermal fluctuations [H El HI El [3, El [9] . 
Our procedure is distinguished by being independent of an interfacial reference plane or 
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presumed symmetry. This feature is particularly useful for studying liquid interfaces near 
non-planar substrates and irregularly shaped solutes like biopolymers. We turn to such an 
application in |2| which shows the instantaneous interfaces of water during a trajectory of 
two hydrated melittin dimers. This example has no obvious spatial symmetries. 

Melittin dimers have hydrophobic domains that are exposed to solvent until they assemble 
to form a stable tetramer. In this assembly, the hydrophobic domains undergo a dewetting 
induced hydrophobic collapse ^TOi [TT| [12] . The rendered interfaces in [2] show how this 
aggregation is highly collective. The concerted motions of water molecules that underlie 
what is pictured would be hard to detect without viewing the instantaneous interfaces. In 
the specific trajectory illustrated, the dimers first come into contact on one end. Collapse 
proceeds through a zipper-like motion during which water is squeezed out of the cavity 
away from the ends of the dimers that have already associated. This process is aided by 
the intermittent formation of vapor-tunnels bridging the unassociated ends of the dimers 
(examples of these vapor tunnels are shown at t = 240ps and t = 377ps in |2]). Due to 
these vapor tunnels, there is an unbalancing of solvent-induced forces, and this unbalancing 
accelerates assembly [131 ttH IIS] • 

Contrasting Mean and Instantaneous Interfaces 

To quantify molecular properties associated with interfaces, it is useful to carry out av- 
erages in terms of the positions and orientations of molecules with respect to the location of 
the interface. In such considerations, the distinction between mean and instantaneous inter- 
faces is significant. This is because interfacial fluctuations can be large, so that a molecule 
located at the position of the mean interface can be often distant from the instantaneous 
interface. Using the mean density or mean interfacial profile to specify the reference surface 
will therefore obscure the true molecular nature of interfacial properties. To illustrate the 
significance, we have examined a few properties associated with the liquid-vapor interface. 
In this case, the mean interface is essentially the fixed reference frame of the Gibbs dividing 
surface [T6] . 

In either frame of reference, we let ai{t) denote the proximity of the ith water molecule 
to the surface. That is, 

ai{t) = {[s(t) - ri{t)] ■ n.{t)}\s(t)=s*it), instantaneous, (5) 
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or 



ai(t) = [(s) - Tiit)] ■ (n), mean. 



(6) 



The angle brackets denote equilibrium average; s*{t) is the point on s{t) nearest rj(t); n{t) 
is the unit vector normal to the instantaneous interface at s(t), i.e., the unit vector in the 
direction of Vp{r,t)\r=s{t)', (n) is the unit vector normal to the mean surface, (s); ai{t) is 
positive if molecule i is on the vapor side of the interface and negative if molecule i is on the 
liquid side of the interface. [While we are here considering specifically liquid-gas coexistence, 
where the mean interface has planar symmetry, the proximity ai{t) defined in Eq. (5) is 
more generally applicable.] 

With this notation, the mean density profiles with respect to either the instantaneous 
water-vapor interface or the mean water-vapor interface is 



the density profile when is in reference to the instantaneous interface, Eq. (5), with that 
when is in reference to the mean interface, Eq. (6). The profile relative to the fluctuating 
instantaneous interface exhibits oscillations indicative of a layering of the atomistic solvent. 
This layering while not without precedent [HI [13, HSl [H] represents a significant departure 
from the more familiar sigmoidal density profile observed relative to the mean dividing 
surface. It is evident from this result that the liquid-vapor phase boundary is indeed well 
described as a sharp surface (with a width of approximately a single molecular diameter) 
dividing the bulk liquid and vapor phases. In this perspective, proposed by Stillinger [3] and 
later treated extensively by Widom [20J and Weeks |4j, the sigmoidal feature seen in the mean 
solvent density profile (|3| is a manifestation of thermal fluctuations in the position of the 
intrinsic surface. An important consequence arises in the context of extended hydrophobic 
interfaces whose density profiles are often drawn to resemble those of a liquid-vapor interface, 
meaning they including a region of depleted solvent density between the hydrophobic surface 
and the bulk liquid. It has been shown that for extended hydrophobic solutes that interact 
with water through even weak dispersive interactions that the solvent density profile is often 
not sigmoidal but exhibits layering |211 [22l [231 1211 [25], [261 [2lj- This observation has been 
used to conclude [211 [231 [27] that water near such hydrophobic surfaces is not like water 




(7) 



where L is the length of the simulation cell parallel to the mean interface. [3] juxtaposes 
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FIG. 3: The mean density profile n(x) (solid line) and the mean orientational property m(x) (dashed 
line) plotted as functions of the distance x relative to either the Gibbs dividing sm'face (top) or 
the instantaneous interface (bottom). The quantity p is the value of bulk liquid density and m is 
the magnitude of the dipole moment of a single water molecule in the simulation. 

near a liquid-vapor interface. |3] shows that such a conclusion is not justified. In particular, 
the mean water density proximal to the instantaneous liquid-vapor interface is layered. The 
layering thus seen in the mean density near hydrophobic surfaces indicates that dispersive 
attractions are sufficient to locate the water- vapor interface on average adjacent to the 
surface. Without those weak attractions the interface would wander. 

Not only does the layering indicate that the instantaneous s neatly divides the vapor and 
liquid phases, the modulation in the solvent density has a significant effect on the orientations 
of water molecules in the vicinity of the water- vapor interface. This orientational structure is 
washed out when the Gibbs or mean surface is used as the reference. Consider, for instance 
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the dipole moment of the ith water molecule, irij, for which the conditional mean of its 
projection parallel to the instantaneous or mean interface is given by 

1 / 



where either 



or 



rrii = (rrij ■ n) |s=s*, instantaneous (9) 



mj = rrij ■ (n), mean. (10) 

These quantities are plotted alongside their respective density profiles in [3j In the frame of 
reference of the instantaneous interface, oscillations of m{x) with respect to x refiects the 
correlations between molecular orientation and proximity to the interface. These correlations 
are not nearly so evident or interpretable in the frame of reference of the mean interface. 

Indeed, orientational structure is sufficiently vivid in the reference frame of the instanta- 
neous interface that we find it informative to also consider quantities like the joint conditional 
distribution 

P{u, u'\x) cx ^ 6(u- cos(^f 6 (u' - cos(^P)) 6 {ai - x)^ , (11) 
where 6^^^ is the angle between the instantaneous surface normal at s* and one of the two 

(2) 

0-H bond vectors of the ith water molecule, and 6- is similarly defined for the other 0-H 
bond vector of that same molecule. This distribution at several values of x, where proximity 
ttj is defined with reference to the instantaneous interface, is shown in |4} 

The panels of that figure show how water molecules near the instantaneous liquid-vapor 
interface adopt orientations consistent with locally favorable hydrogen bond patterns. In 
the first layer of liquid, where n{x) has its first peak {x ~ 1-7A), the figure shows that 
water molecules most likely adopt orientations where both OH bonds are aligned to donate 
hydrogen bonds to other water molecules residing in that first layer. The geometry of a 
water molecule makes it impossible for a molecule in the first layer to donate in two such 
hydrogen bonds when x < O.SAand x > 2.6A. Due to this constraint, the population of 
water molecules at a; ~ O.SAmost likely orient only one 0-H to donate hydrogen bond to the 
first layer of liquid while the other OH bond pointed into the vapor phase {£ ^ 1). Similarly, 
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FIG. 4: The joint probability distribution P(u,u'\x) governing the orientations of each OH bond 
of water molecules with x = 0(a), x = 0.8A(b), x = 1.2A(c), and x = 2.oA(d). A value of i = 1, 
I = 0, oi i = —1 refer to OH bonds that are oriented parallel, perpendicular, or anti-parallel to 
the local surface normal respectively. The normalizing factor P(u,u' \oo) is the value of the joint 
probability distribution in an isotropic environment. 
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water molecules near x ~ 2.6Amost likely donate one hydrogen bond into the first layer and 
another into the second layer. 

While not illustrated explicitly here, we find that orientational structure extends from the 
instantaneous interface into the bulk liquid as far as SA. For x > SA, the solvent orientational 
structure ceases to exhibit the infiuence of the interface. 



Methods 



The numerical simulation study consisted of SPC/E water molecules [28] in a 36 x 36 x 
lOOA'^ simulation cell periodically replicated in the x and y Cartesian directions. At a tem- 
perature of 298 K the 1387 water molecules in the simulation cell form a liquid slab spanning 
the periodic boundaries that is approximately 36A thick (in the z Cartesian direction). The 
liquid-vapor phase boundaries present in the simulation serve as a natural barostat and 
therefore the liquid can be regarded as a system being held at constant pressure. Electro- 
static interactions are treated with two-dimensional Ewald summation [29], and molecular 
constraints are enforced with the RATTLE algorithm [30]. Statistics were generated through 
six independently equilibrated molecular dynamics simulations each run with a time step of 
2.4fs for about Ins with nuclear coordinates written out every 50 time steps. 

To identify the water-vapor interface the density field p(r, t) was computed on a cubic 
lattice with a lattice spacing of lA. For the spatial coarse-graining of the density field, 0(r; S,) 
was truncated and shifted to be both continuous and zero at a distance of 3^. The Gaussian 
width ^ was selected by first computing a measure of the average amount of interfacial area 
in the system, 

/ 0[P(r,^)-c]e[c-p(r + /z,t)]dr^, (12) 
where Q{x) is the Heaviside function which is given by. 



Q{x) 



if a; > 0, 
if X < 0. 



We have found that a value of / = lAis sufficiently small to ensure an accurate measurement 
of A. For small values of ^ the quantity A decreases with increasing ^, eventually reaching 
a constant value oi A = the projected area of a single planar interface. For large values 
of ^, A = L^, the projected area of a single planar interface. For smaller values of ^, 
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A> L"^, which happens when the density field p{r,t) in the hquid contains cavities and/or 
the planar interface develops overhangs. We find that a value of ^ = 2.4Ais large enough 
to essentially eliminate the occurrence of interface overhangs and bubbles within the liquid 

phase. Computing the entire density field in this manner for a single configuration (time 
step) took 3.4 seconds on a single processor of a modern desktop computer. 
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